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The finite amplitude metliod (FAM), wliicli we have recently proposed (T. Nakatsukasa, T. In- 
akura, and K. Yabana, Phys. Rev. C76, 024318 (2007) ), simplifies significantly the fully self- 
consistent RPA calculation. Employing the FAM, we are conducting systematic, fully self-consistent 
response calculations for a wide mass region. This paper is intended to present a computational 
scheme to be used in the systematic investigation and to show the performance of the FAM for 
a realistic Skyrme energy functional. We implemented the method in the mixed representation 
in which the forward and backward RPA amplitudes are represented by indices of single-particle 
orbitals for occupied states and the spatial grid points for unoccupied states. We solve the linear 
response equation for a given frequency. The equation is a linear algebraic problem with a sparse 
non-hermitian matrix, which is solved with an iterative method. We show results of the dipole 
response for selected spherical and deformed nuclei. The peak energies of the giant dipole resonance 
agree well with measurements for heavy nuclei, while they are systematically underestimated for 
light nuclei. We also discuss the width of the giant dipole resonance in the fully self-consistent RPA 
calculation. 

PACS numbers: 



I. INTRODUCTION 

Nuclear density-functional theory has been successful 
for a systematic description of ground-state properties 
such as binding energies, density distributions, deforma- 
tions, and so on. In the last decade, systematic density- 
functional calculations for a whole nuclear chart have 
come to be an ordinary work where a variety of energy 
functionals are employed One of the major goals 

of such systematic investigations is an assessment and an 
improvement of the energy functionals. A rapid progress 
of experimental research on unstable nuclei also necessi- 
tates a development of an accurate and universal energy 
functional that describes nuclei far from the stability line. 
These efforts ultimately aim to establish a theory with a 
high predictive power and high accuracy. 

The nuclear density-functional theory is also capable 
of describing excited state properties. The random-phase 
approximation (RPA) , which is derived by linearizing the 
time-dependent mean-field equation , describes the nu- 
clear excitation as a small-amplitude oscillations around 
the ground state. It has been successful for describ- 
ing both giant resonances and low-lying modes of ex- 
citations. In practical calculations, the RPA equation 
has been solved mostly in the matrix diagonalization 
scheme. In these calculations, it was common to ig- 
nore a part of the residual interaction, sacrificing the 
full self-consistency. This is because the inclusion of the 
full residual interaction for a realistic interaction involves 
a cumbersome and complicated task. Recently several 
groups have reported fully self-consistent RPA calcula- 
tions [1, H, H, 0, [1| . However, they are mostly restricted 
to spherical nuclei. There are a few recent attempts for 

deformed nucle i d [Mini, [3 in, El- 

Recently, we have proposed a new method to solve the 



RPA equation, the finite amplitude method (FAM) [T^. 
The FAM allows us to evaluate the fully self-consistent 
residual fields as a finite difference, employing a compu- 
tational code for the static mean-field Hamiltonian alone 
with a minor modification. In Ref . [l5| . we showed that 
the FAM works accurately in solving the RPA equation 
for a deformed nucleus ^''Ne with a simplified interaction. 

Employing the FAM, we are undertaking a systematic 
investigation of nuclear response for a wide mass region 
{A < 100). The purpose of this paper is to explain the 
implementation of the FAM with a realistic Skyrme func- 
tional to be employed in such systematic investigation. 

Our scheme is summarized as follows: 

1. Fully self-consistent RPA calculation. The iden- 
tical Skyrme energy functional is used for both 
ground state and linear response calculations, in- 
cluding spin-orbit, Coulomb, and time-odd-density 
terms. 

2. Applicable to deformed nuclei. To describe any 
kind of deformation, we employ a three-dimensional 
(3D) real-space grid representation. 

3. Linear response equation is solved employing an 
iterative method for the external field with a given 
frequency. 

We employ a mixed representation in which the RPA 
forward and backward amplitudes are described with or- 
bital indices for occupied single-particle states and the 3D 
grid-points for unoccupied states. In the mixed represen- 
tation, the RPA equation is a linear algebraic equation 
with a sparse, non-hermitian matrix. The dimension- 
ality is typically an order of 10®. We will examine and 
compare performance of different solvers for the linear al- 
gebraic problem, and demonstrate that the FAM works 
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excellently for realistic calculations with complicated en- 
ergy functionals. 

At present, we ignore the pairing terms. The pairing ef- 
fects will not be important for high-lying negative-parity 
excitations such as the dipole excitation which we investi- 
gate in this paper. However, the pairing may have signif- 
icant effects on the low-lying modes such as quadrupole 
excitations. An extension to include the pairing is now 
under study. In order to treat the particle escape , w e 
need to impose a continuum boundary condition 
Although it can be managed in the present scheme 
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the accurate description requires a heavy computational 
task. For a systematic investigation which we pursue, 
we must find a compromise between the accuracy of the 
calculation and the computational feasibility. We will ap- 
proximately take into account the continuum effects by 
solving the RPA equation inside a 3D spherical box of a 
large radius. 

This paper is organized as follows. In Sec. [TTl we re- 
view formalism of the FAM. The numerical details and 
the related techniques are discussed in Sec. Ill CI To 
show the accuracy of the method, the FAM results are 
compared with the fully self-consistent RPA results in 
the conventional diagonalization scheme. Calculated re- 
sults are shown in Sec. IIIII for selected nuclei in light 
and heavy mass regions. Convergence with respect to 
the model space, dependence on the functionals will be 
discussed. Finally, our fully self-consistent RPA results 
are compared with experimental data. Conclusions are 
presented in Sec. |IVl 



II. LINEAR-RESPONSE CALCULATION WITH 
A SKYRME FUNCTIONAL 

A. Linear response equation in the mixed 
representation 

Under a weak, time-dependent external field Vcxt{t), 
the transition density Sp{t) follows the equation (h ~ I) 

i^^Spit) = [ho , dp{t)] + [Kxt(0 + Sh{t) po] , (1) 

where po and /iq — h[pQ] are the ground-state den- 
sity and the single-particle Hamiltonian, respectively. 
dh(t) is a residual field induced by density fiuctuation, 
h[p{t)] — ho + 5h{t). Assuming that 5p{t)^ Kxt(i), 
and 5h{t) oscillate with a frequency lj like 5p{t) = 
(5p(w)e~^'^* -I- (5p'f(w)e*"*, Equation ^ is recast to 



uj5p{uj) = [ho, 5p{uj)] + [Vc^t{i^) -f- 5h{uj),po] 



(2) 



Since 5p{uj) is not necessarily hermitian, we need for- 
ward and backward amplitudes, [Xiiuj)) and [Yi[uj)), to 
express the transition density 6p{uj). 



(3) 



where A = N + Z is the mass number of the nucleus and 
the orbitals 4>i are occupied orbitals in the ground state, 
hol'f'i) — ^il'Pi) (* = Ij"'' i^)- Substituting this ex- 
pression into Eq. ^ , we obtain the RPA linear-response 
equation: 

uj\Xi{uj)) = {ho - a) \X^{uj)) 

+P{V,M+Sh{Lo)}[(b,), (4) 
-u{Y,{uj)\ = {Y,iio)[iho-e,) 

+ (0»|{K=xtM+<5/iH}P, (5) 

where P denotes the projector onto the particles space, 

P = 1 - Ef=i \(t>^){<P^\■ In Eqs. O and (P, if we ex- 
pand 5h(uj) with respect to the forward and backward 
amplitudes in the linear order, we obtain the well-known 
A — B matrix form of the RPA equations [1, [lB| ■ When 
the single-particle Hamiltonian and the external field are 
both local in coordinate space, we may write Eqs. ^ and 
([5]) conveniently in the coordinate representation, 

+ P{V,^titu;) + 6hitio)}MO, (6) 
-u;Y:{^,u;) = po(e)-eOr,(C,c.) 

+ P {vUt ^) + 5h\t u)] * 17) 

where the coordinate ^ may contain the spin index a as 
well as the spatial coordinate r,£^ = (r, a). This is of- 
ten referred to the mixed representation [l^, [H, [13 • It 
should be noted that, since Sh{uj) has a linear depen- 
dence on Xi{oj) and Y*{(jj), this is an inhomogeneous 
linear algebraic equation of the form Ax = b where 
X = {Xi{£_,uj),Y* {£^,uj)) [l^. In our implementation, 
we employ the grid representation of the 3D Cartesian- 
coordinate space. Denoting the number of the spatial 
grid points as Np, the dimension of the vector x is 
Np X A X 2 X 2 (a factor of two for the forward and 
backward amplitudes and another factor of two for the 
spin degrees of freedom). Nf is order of 10,000 (see Sec. 
Ill C II) . For systems described by local potentials, we 
need not introduce any further truncations in the particle 
space. In the 3D grid representation the treatment of the 
particle escape is better controlled than that in other rep- 
resentations, for example, the harmonic-oscillator-basis 
representation. 



B. Finite amplitude method (FAM) 

If we write the residual field 6h{£^,u!) explicitly with 
the forward Xi{^,uj) and backward Yil^jUj) amplitudes, 
it requires us to calculate the residual two-body kernel, 
v{^,^') = Sh{^,uj)/6p{£,'). However, if one employs a re- 
alistic energy functional, the explicit construction of the 
residual two-body kernel involves a complicated coding 
and a large numerical task. For this reason, it has been a 
common procedure to approximate or ignore a part of the 
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residual interaction in most RPA calculations. For exam- 
ple, the velocity-dependent terms are sometimes replaced 
with the Landau-Migdal interaction. The spin-spin, the 
spin-orbit, and the Coulomb residual interactions are of- 
ten neglected. 

To achieve a fully self-consistent RPA calculation, we 
believe that the FAM [l5| is one of the simplest methods, 
at least, with respect to programming. In the FAM, we 
evaluate Sh{£^,uj) directly from the ground state Hamil- 
tonian h[p\ by the finite difference method, thus avoiding 
an explicit construction of the residual two-body kernel. 
In this section, we recapitulate the essence of the FAM. 
Readers are referred to Ref. [l^ for details. 

The residual field 6h{£^,uj) depends linearly on 5p{uo): 
5h{ijj) — dh/ dp\p=pa ■Sp{uj). In the first oder with respect 
to a small parameter rj, we have 



ho + ri5h[uj) ^h[po + ri5p{uj)\ = h 



(8) 



Here, the one-body pseudo-density operator Pni^) = Po + 
ri6p{u>) is a non-hermitian operator expressed with bra 
and ket vectors as follows; 

p^iu;) = d'^') + ^I^^M)) (('^^l + v{Y^{u;)\) . (9) 



In the linear order in ?/, Equation ([8]) is expressed as 



dh{uj) = -{h [pr,{u})] - ho) . 
V 



(10) 



This is the basic result of the FAM. Equation (fTU)) indi- 
cates that h [p^{uj)] can be evaluated in the same manner 
as in ho = h[po], with a replacement of the ordinary den- 
sities with the pseudo-densities: For instance, the pseudo- 
local- density, pseudo-kinetic-density, and pseudo-current- 
density are given by 

i a 
i a 



(11) 



respectively, where '^i — (^i ^ r\Xi(ijj) and -0^ = 0^ + 
r]Yi(ijj). The spin-dependent pseudo-densities are also de- 
fined in the same manner. Once Xiiio), and r/ are 
given, the pseudo-densities are calculated as in Eq. pTjl . 

Equations ([6]) and ([7]) are the linear algebraic equation 
of the form, AS" = h. Here, 



-Fext(e,^)0^(e) 



(12) 

To solve Eqs. dH) and ©, we utilize the iterative scheme, 
as described below. In the iterative scheme, we need not 



construct the matrix elements of A explicitly but always 
evaluate Ax employing the FAM for a given vector x. 



Ax = 



{{hoiO - e,+u;*)m,uj) + Sh^{^,uj)MOy 



(13) 

where Sh{£^,uj)(j)i{^) is calculated using Eq. (fTO|) and 
(5/i^(^, a;)(/)j(^) using the same equation but replaced 
by pt. To find a solution of Eqs. ([6]) and (O, we em- 
ploy extensions of the conjugate gradient method for lin- 
ear algebraic equations involving a non-hermitian matrix. 
In the literature, quite a few solvers for linear algebraic 
equations involving a non-hermitian matrix have been 
developed. However, we find only a few methods work 
for the present problem. We will discuss it in Sec. Ill Cl 

In order to calculate the strength function for a given 
one-body operator F, we adopt an external field of 
Vc^t{t) = rjFe-"^* -\- r]*Fh"^K Then, the transition 
strength is expressed with the forward and backward am- 
plitudes. 



^^W^ ^ -^Im^{(0,|Ft|X,(a;)) + (r,(c.)li^t|^^)} 

i 

= J2\{n\F\0)\^SiE-E^), (14) 

n 

for a real frequency uj = E. Here, \n) are energy eigen- 
states of the total system. In the following calculations, 
we use complex frequencies with a finite imaginary part, 
Lo = E -\- 17/2. Then, the transition strength becomes 



dB{E;F) 
dE 



7 |(n|F|0)p 
2n^\iE-E.,f + i^//2f 

|HFt|0)p 1 



{E + E^y + {j/2y 



(15) 



The second term in the right hand side is not important, 
if 7 <C i? + En- For an hermitian operator F, this leads 
to 



dB{E; F) _ 2E'y 
dE ~ ~V 



E 



En\{n\F[ 



E^ - El 



^272 



(16) 



where Ef^ = El+ 72/4. 

In this article, we consider an electric dipole operator 
for F. 



N ^ Z ^ 

p— 1 n— 1 



(17) 



and similar operators for D^^ and Dy^ . The photoab- 
sorption cross section in the dipole approximation is 
given as follows [l8| , 
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C. Numerical details 

1. Adaptive 3D grid representation 

We employ a model space of 3D grid points inside a 
sphere of a radius i?box- AH the single-particle wave func- 
tions and potentials except for the Coulomb potential are 
assumed to vanish outside the sphere. For the calculation 
of the Coulomb potential, we follow the prescription in 
Ref. [19]. The differentiation is approximated by a finite 
difference with the nine-point formula. 

In order to obtain nuclear response with a reasonable 
accuracy, we need a large model space, typically i?box ^ 
15 fm (see Sec. IIIIB This is much larger than the 
typical nuclear size Rq- Outside the nuclear radius Ro, 
we do not need to use a fine mesh. Thus we adopt an 
adaptive coordinate system to reduce the number of grid 
points in the outer region, Rq < r < i?box- We use the 
following coordinate transformation, (w, u, w) — > (x, y, z), 
which was also adopted in Ref. , 



x{u) — ku 



sinh(w/xo) 



(19) 



and the same form for y{v) and z{'w). In this transforma- 
tion, we have x{u) « u for spatial region of u <gc xq, and 
x{u) « ku for xq. We adopt a uniform mesh spacing 
of A/i — 0.8 fm in the {u,v,w) space. The parameters, 
k — 5.0, Xq — 8 fm, and n — 2, are used in the following 
calculations. The number of grid points in the sphere of 
r < Rhox is significantly reduced; 27,609 — > 11,777 for 
-Rbox = 15 fm and 221, 119 39, 321 for i?box = 30 fm. 



2. Choice of Iterative algorithms 

In this section, we discuss performance of different it- 
erative solvers to solve the linear response equations. We 
first discuss the calculation for the ground state. To ob- 
tain ground-state solutions, the imaginary-time iterative 
method is used We impose constraints in the itera- 
tion process so that the center-of-mass always coincides 
with the origin of the coordinate system and the princi- 
pal axes of the density distribution coincide with three 
Cartesian axes of x, y, and z. It is important to obtain 
a well converged ground state solution, since the conver- 
gence properties of the linear-response calculation cru- 
cially depends on this. To assure the strict convergence 
in the ground state calculation, we impose the conver- 
gence condition not only for the total energy, but also 
for the deformation parameters (/3,7), single-particle en- 
ergies Ci, single-particle angular momenta, and so on. 

For the Skyrme energy functional, the RPA matrix A 
which appears in Eqs. ^ and ([7]) is sparse in the r-space 
grid representation. Therefore, the iterative methods, 
such as the conjugate gradient (CG) method [IJ], should 
work efficiently. The CG method is very powerful for the 
hermitian matrix. However, since we calculate for the 
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FIG. 1: Convergence property of different iterative methods 
to solve the linear response equations (|6]) and ^ for electric 
dipole response in ^^O calculated at complex frequencies of 
= -I- 0.5i MeV (left) and 10 + 0.5i MeV (right). Relative 
residue, r„, is shown as a function of iteration number n. 



complex frequency lu, the matrix A is no longer hermi- 
tian. Therefore, a simple CG method is not applicable. 
There are a lot of variants of the CG method extended 
for non- hermitian problems. We test some of them: Bi- 
conjugate gradient (Bi-CG) method [2l|, generalized con- 
jugate residual (GCR) method [25|, generalized product- 
type bi-conjugate grad ient (GPBi-CG) method [1^, Bi- 
CGSTAB method^, and Bi-CGSTAB2 method H]. 
In the original GCR method, it is necessary to store all 
the vectors xq, - ■ ■ ,2^n-i m order to calculate the vec- 
tor at the n-th iteration. However, this is impractical 
because it requires heavy computational task and huge 
memory size. To avoid the problem, we restart the GCR 
procedure every twenty iterations. In Fig. [U we show 
the convergence behavior of the different iterative solvers. 
The magnitude of the relative residue, 



r„ = \b- Aa;7,|/|6| 



(20) 



is plotted against the number of iterations. 

At very low energies {uj = + 0.5i MeV), all the 
solvers except for the Bi-CG method provide quickly con- 
verged results. On the other hand, at higher energies 
(w = 10 + 0.5i MeV), only the GCR and the GPBi-CG 
methods reach the convergence. In most cases, the con- 
vergence of the GPBi-CG method is faster than the GCR. 
However, after the convergence, the residue by the GPBi- 
CG method start to increase again. Therefore, only the 
result of the GCR method shows a stable convergence 
property in Fig. [1] We thus employ the GCR in the fol- 
lowing calculations, although it has a disadvantage that 
it requires larger computer memory size than other meth- 
ods based on the Bi-CG. 

In the iteration procedure, we need to set up the ini- 
tial vector a'o • It turns out that the convergence property 
depends very little on the selection of the initial vector. 
We simply take the initial vector of xq = 0. We stop 
the GCR iteration procedure if the relative residue be- 
comes smaller than a threshold, r„ < 10~^. We find 
the typical number of the iteration necessary to reach 
the convergence is about lO'^ around the excitation en- 
ergy region of giant dipole resonance (GDR), while the 
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convergence is much faster for smaller E{= Rew). The 
iteration number to reach the convergence also depends 
on the imaginary part of the frequency, 7(= 2Imci;), and 
is roughly proportional to I/7. 

We will later show how the convergence depends on 
other parameters. In the inset of Fig. [3l we show the 
residue, r„ of Eq. as a function of the CPU time 

for calculations employing different spatial size i?box- We 
can see that the computational time scales linearly both 
with the number of grid points and with the particle num- 
ber A. 

In Sec. Illlj the calculation is performed with the fol- 
lowing settings, unless otherwise specified: The energy 
range of < E < 38.1 MeV, discretized in spacing of 
AE = 0.3 MeV (128 points). The imaginary part of the 
frequency is fixed at 0.5 MeV, corresponding to 7 = 1.0 
MeV. The calculated strength is interpolated using the 
cubic spline function. In most cases, we have utilized 
PC cluster systems with either 64 or 128 processors in 
parallel. 

3. Choice of the FAM parameter rj 
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E [MeV] 

FIG. 2: Comparison of the results for electric dipole response 
of ^^O calculated by the present linear response formalism 
with FAM (dots) and by the eigenvalue formalism with an 
explicit construction of the residual interaction (curve) . Both 
calculations use the cubic box of (21 fm)"^ and the SIII func- 
tional. See text for details. The vertical lines indicate eigenen- 
ergies and strengths of the RPA normal modes calculated by 
the matrix diagonalization method, whose vertical magnitude 
is in units of mb-MeV. 



The FAM parameter 77 in Eq. (fTO| should be as small 
as possible to validate the linearity. In practice, there is a 
lower limit to avoid the round-off error. In the numerical 
computation, we use real variables with double precision 
(8 bytes) and complex variables with 8x2 bytes. We use 
the following value for p which differ for every iteration 
to ensure the linearity [l5[: 



max{iV(A:),iV(V)}' 



iV(^) 



1 

A' 



(21) 

with 77Lin = 10 The convergence property of the GCR 
iteration is not sensitive to the value of rjhim as far as the 
value is in the range of ?7Lin = 10^^ — 10^^. However, if 
we use too small value, r„ oscillates in the iteration and 
never reaches the convergence. 



III. CALCULATED RESULTS 
A. Numerical accuracy 

In this section, we demonstrate that the FAM, de- 
scribed in Sec. [TTl really works as a simple and accurate 
method for the fully self-consistent linear-response calcu- 
lation with realistic Skyrme functionals. In Ref. [lH|, we 
have already done it for a simple Bonche-Koonin-Negele 
interaction [2q . 

In the following calculations, we use the Skyrme en- 
ergy functional. We adopt the same functional form as 
that of Ref. [l^j (Appendix A). In Ref. [l^l, every single- 
particle orbital is assumed to have a definite parity and 
z-signature symmetries. However, we do not assume any 



of these in this paper. The functional depends on den- 
sities Prir), kinetic densities Tr{f), spin-orbit densities 
Jri^), current densities jri^), and spin densities Prir), 
where t — n and p. The exchange part of the Coulomb 
energy is evaluated by the Slater approximation. The 
center-of-mass correction is achieved by using the recoil 
nucleon mass, M — > AM /{A — 1). 

In Fig. [21 we show photoabsorption cross section for 
^^O calculated by two different methods: The dots are 
obtained by solving the linear response equations ^ and 
^ in which the FAM is employed to evaluate the resid- 
ual field. The vertical lines are obtained by the matrix 
diagonalization method [l^, in which the residual in- 
teraction taking a second derivative of the Skyrme en- 
ergy function is constructed explicitly. The positions of 
the vertical lines indicate the energy eigenvalues, while 
the heights show the calculated cross section. The solid 
curve is obtained by smoothing the vertical lines with the 
energy-weighted Lorentzian. In both calculations of the 
linear response and the matrix diagonalization, the same 
Skyrme functional and the same grid representation of 
the r-space are used. The width parameter is set com- 
mon, 7 = 2 MeV. In the matrix diagonalization method, 
the spurious solutions corresponding to the center-of- 
mass motion appear around a few tens keV. They are 
not shown in the figure. 

This figure clearly demonstrates that the linear re- 
sponse calculation with the FAM denoted by dots com- 
pletely agrees with the fully self-consistent RPA calcula- 
tion with an explicit construction of the residual inter- 
action. We should note that the computer program of 
the diagonalization method requires heavy and compli- 
cated coding for the residual interaction, while the pro- 
gram coding with the FAM can be achieved with the 
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FIG. 3: Photoabsorption cross sections in ^''O and '*''Ca cal- 
culated with the FAM in a spherical box of different sizes: 
i?box = 10 fm, 15 fm, 20 fm, and 30 fm. "ABC" indicates 
a result with the absorbing boundary condition which prop- 
erly treats the continuum effect [Til ]. The SkM* interaction 
and the smoothing parameter 7=1 MeV are used. Inset: 
Convergence properties as a function of the CPU time. 



static mean-field calculation alone, with a minor modifi- 
cation. Thus the accuracy and the simplicity of the FAM 
is clearly demonstrated. 

Let us next consider the merit of the linear response 
calculation for a fixed frequency in comparison with the 
matrix diagonalization method, apart from the FAM. In 
the matrix diagonalization method, the number of nor- 
mal modes increases rapidly as the excitation energy be- 
comes higher. Even when we assume the parity and the 
^-signature symmetries, the number of excited states is 
an order of one thousand in the energy range below 50 
MeV. If we achieve matrix diagonalization in the mixed 
representation, we must calculate eigenvalues and eigen- 
vectors of the RPA solution one by one. Then the cal- 
culations of a thousand of eigenstates require huge com- 
putational time. Therefore, to calculate responses for a 
wide energy region, the computation time with the linear 
response method is much smaller than that of the diago- 
nalization method. Furthermore, we can easily parallelize 
the calculation in the linear response calculation, assign- 
ing different processors to solve the equation at different 
frequencies uj. 

Combining these features, we consider the linear re- 
sponse calculation with the FAM is the best choice for 
the systematic, fully self-consistent calculations over a 
periodic table. 



B. Light nuclei 

1. Dependence on box size and smoothing parameter 

We show how the results depend on the box size i?box 
and the smoothing parameter 7. In the present mixed 
representation, the model space is specified with i?box 
and the 3D mesh size Ah. The dependence on the mesh 
size is discussed in details in Ref. 
low-energy excitations in ^^O and 
this study. Ah « 0.8 fm produces a converged result 
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FIG. 4: Photoabsorption cross sections in spherical nucleus 
^''Ca calculated with different Rhox and with different 7 val- 
ues. The functional parameter set of SkM* is used. 




12| for spurious and 
^°^Pb. According to 



FIG. 5: Photoabsorption cross sections in deformed nucleus 
^*Mg calculated with a spherical box of iibox = 10 fm, 15 
fm, 20 fm, and 30 fm. Dashed, long dashed, and solid lines 
correspond to cross sections of iC" = mode, K = 1 mode, 
and the total one, respectively. The functional parameter set 
of SkM* and the smoothing parameter 7 = 1 MeV are used. 



In addition, we have also carefully examined the conver- 
gence with respect to the grid spacing in the adaptive- 
coordinate. We adopt xq — 8 in Eq. p^ . 

In Fig. [3l results of the calculations with i?box — 10, 
15, 20, 25, and 30 fm are shown. As is seen in Fig. [31 
small peaks in higher energies {E > 20 MeV) are dis- 
appearing with increasing i?box- These small peaks at 
high energies are all spurious, due to the discretization 
of the continuum. In fact, the calculation with the proper 
continuum boundary condition p^ . which is denoted as 
"ABC", shows only a smooth tail in this energy region 
above the main GDR peak. With the energy resolution of 
7=1 MeV, the effect of the discretized continuum is still 
visible at i? > 20 MeV, even for the case of i?box — 30 
fm. In order to remove the spurious effect, we need in- 



TABLE I: _Rbox-dependence of the GDR peak positions and 
widths in units of MeV, for ^''O and "^^Ca. The smoothing 
parameter of 7 = 1 MeV is used. See text for details. 
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-^box 


= 10 fm 


18.993 
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18.798 


3.876 


-^box 


= 15 fm 
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-^box 


= 20 fm 


19.141 
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19.151 


3.435 




= 30 fm 
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3.035 


19.165 


3.330 






dB{E)/dE 


0"abs(i5) 


^"Ca 


-f/pcak 
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-fi^pcak 


r 


^box 


= 10 fm 


17.794 


1.545 


17.766 


1.550 


Rhox 


= 15 fm 


17.797 


2.334 


17.796 


2.492 


-^box 


= 20 fm 


17.904 


2.789 


17.901 


2.886 


-^box 


= 30 fm 


17.963 


2.956 


17.959 


3.014 



crease either i?box or 7. In Fig. [H we show results for 
'^'^Ca calculated with different box sizes and with different 
magnitudes of the smoothing parameter 7. Using 7 > 2 
MeV, the calculation roughly converges with i?box > 15 
fm. Thus, in the calculation with the vanishing (box) 
boundary condition, it is difficult to distinguish small 
physical peaks from unphysical ones in the high-energy 
continuum region, E > 20 MeV. In this article, we con- 
centrate our discussion on gross structure of the main 
GDR peaks. 

Next, in Fig. O we show results for ^^Mg. Since the 
ground state of ^^Mg is deformed with a prolate shape 
{P2 — 0.39), the GDR splits into two major peaks, K — 
mode (dashed line) and K = 1 mode (long dashed line). 
The K — mode corresponds to the oscillation parallel to 
the symmetry axis (z direction) and K = 1 mode to that 
perpendicular to the symmetry axis (x and y directions). 
In the calculations with i?box — 10 fm and 15 fm, we 
observe a peak at E = 25 MeV for the K = 1 mode. 
However, this peak becomes smaller for larger i?box- This 
indicates that the peak at 25 MeV comes from the effect 
of the box discretization and is spurious. 

Although the calculated profile functions of cross sec- 
tion in the high-energy continuum region are affected by 
the box discretization, the gross property of the GDR 
is less sensitive to the value of i?box- We estimate the 
GDR peak energy and the width according to the electric 
dipole {El) strength distribution and the photoabsorp- 
tion cross section (oscillator strength distribution). For 
spherical nuclei, the nuclear response does not depend on 
the direction of the external dipole field. Thus, we can 
arbitrary choose it, for instance, F = Df^. Then, El 
strength is fitted by a single Lorentzian function. 



dB{E-D^^) 



dE 



L{E) = \ 



i Btot X L{E), 

r/2 



-^{E^ iJpeak)' + (r/2)2 

where -Epcak and T are fitting parameters, and i?tot 



(22) 
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FIG. 6: Estimated GDR width for '*'^Ca as a function of the 
parameter 7. 



J dEdB{E)/dE. The fitting is also performed for the 
photoabsorption cross section <t{E), using the energy- 
weighted Lorentzian, ftotxEL{E). The results are shown 
in Table U These two definitions of the peak energy and 
the width lead to nearly identical results. 

The dependence of -Epcak on the box size i?box is very 
small. On the other hand, the GDR width F shows a 
substantial dependence on the box size. For 7=1 MeV, 
to estimate the width with a reasonable accuracy requires 
^box > 20 fm. However, the width also depends on the 
smoothing parameter 7, which is shown in Fig. [6l The 
value of F is converged at i?box > 20 fm for 7 = 1 MeV. 
When we use a larger 7 value, the model space of i?box = 
15 fm is good enough to calculate the width. From the 
approximate linear dependence of F on the parameter 7, 
we can estimate the width at 7 = as F^^^ sa 2 MeV. 
This value can be regarded as the damping width in the 
RPA, namely due to the Landau damping and escaping: 
pRPA ^ pLandau _|_ pT^ Siucc the Smoothing parameter 
7 can be regarded as the spreading width, 7 = F^, we 
may extract from experimental data. This will be 
discussed in the foUowings. 



2. Comparison with experiments and Skyrme-functional 
dependence 

Calculated photoabsorption cross sections for spheri- 
cal nuclei, ^^O, ^"Ca and deformed nucleus '^^Mg, are 
compared with experimental data in Fig. [T] The calcu- 
lated energies of the GDR peaks are underestimated by a 
few MeV, but the overall structures are reproduced. For 
spherical nuclei, the GDR widths calculated with 7 = 1 
MeV are narrower than the corresponding experimental 
data. This seems to suggest that the spreading width 
F-^, which takes account of effects decaying into com- 
pound states, such as two-particle- two-hole excitations, 
is slightly larger than 7 = 1 MeV. The cross section at 
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FIG. 7: Calculated (upper panels) and experimental (lower 
panels) [2^, [2^ photoabsorption cross sections in spherical 
nuclei '''O, ■'"Ca and deformed nucleus ^*Mg. We use the box 
size -Rbox = 30 fm and the SkM* parameter set with 7 = 1 
MeV. 
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FIG. 9: Left panels: Contour maps of the calculated transi- 
tion densities oi K = Q mode, Im5p(a;), in the y-z {x = 0) 
plane aX uj — 19.4 + 0.5i MeV in ^^O (upper panels) and 
16.1-|-0.5j MeV in ^^Mg (lower). Contour lines with positive, 
negative, and zero values are denoted by solid, dashed, and 
long-dashed lines, respectively. Thick solid lines correspond 
to nuclear surface, defined by / {y^) + j {7?) — 1. Right 
panels: Transition densities predicted by the Tassie model. 
The magnitude is fitted to the RPA result at the maximum 
value. See text for details. 



FIG. 8: Comparison of the calculated photoabsorption cross 
sections with different Skyrme parameter sets and experimen- 
tal data in and ^*Mg. The calculation is performed with 
the box size .Rbox = 25 fm. 



higher energies is larger in experiment. This may be due 
to effect of the ground-state correlation produced by the 
tensor and short-range parts in the two-body nuclear in- 
teraction [13. For ^^O, all the calculations fail to re- 
produce a double-peak structure of the GDR, present in 
experimental data. If we neglect the time-odd spin densi- 
ties, the double peaks can be reproduced [13] ■ The GDR 
peak energy in ^"Ca is better reproduced compared to the 
case for ^^O. The calculated cross sections in the high- 
energy region also becomes closer to the experimental 
data. It seems that the discrepancy is more prominent 
for lighter nuclei. 

For the deformed nucleus ^''Mg, the GDR peak split- 
ting caused by the ground-state deformation well agree 
with the experiments, although the magnitude of the de- 
formation splitting is slightly too large in the calcula- 
tion. We may interpret that the experimental GDR peak 
around i? = 20 MeV is associated with the -fC = mode, 
and those at -E — 22 ~ 25 MeV correspond to the K = \ 
mode. The double-peak structure of the K = \ GDR 
peak is well reproduced. Approximately, the calculated 
cross section is shifted to lower energy from the experi- 
mental ones, by about 3 MeV. 

We next show results obtained with different Skyrme- 



parameter sets in Fig. % SkM* [l^l, SIII [13, SGII [H, 

and SLy4 Two parameter sets, SkM* and SLy4, 

produce similar results. In ^^O, the GDR peak posi- 
tion is calculated around -E = 19 MeV. In contrast, the 
SIII functional produces the GDR peak near E = 21 
MeV. The result with the SGII functional is interme- 
diate, around E — 20 MeV. However, all the calcula- 
tions underestimate the experimental GDR peak energy, 
E 22 MeV. This discrepancy on the GDR peak en- 
ergy will disappear for heavier nuclei (see Sec. IIII C[) . 
The cross section at higher energies E > 25 MeV is also 
systematically underestimated in the calculations. For 
deformed ^"^Mg, all the calculations again underestimate 
the peak energies. Both K = and K = 1 peaks 
obtained with the SIII functional are located highest in 
energy compared to the others. 

3. Transition density 

We obtain the forward and backward amplitudes, 
Xi{^,uj) and Fi(^,a;), using the FAM. Then, the lo- 
cal part of the transition density in the intrinsic frame, 
Eq. ((31), can be calculated as 

Sprir,u;) = {0*(C)X,(e,c.) + Y*{tu;)MO} , 

(23) 

where r = n or p. In this section, we adopt the external 
field of L»f 1 {K = 0). Using Sp{uj) with uj = E + ij/2, 
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the strength function can be written as 



dE TT 



2E-f 



E 



Ze 
~A 



zSpn{r,uj] 

2 



(24) 



Substituting {n\D^ \0) by the following expression 

(n|Df 10} = 1 dr^z{n\pp{m ^z{n\p^{r^)\0), 

(25) 

we find 

. . , ^ K(0|i^f |n)(n|Pn(rt(f)|0) 
Im dp„(p) (r, oc 2^ ; —^2 ' (26) 



E^ - El 



where Prir) — X^iGT*^!^^ ^0- If there is a single state 
\n) that gives a dominant contribution in the vicinity of 
E = En (within the width of 7), we have Im Sp{uj) oc 

The left panels of Fig. [5] show the calculated tran- 
sition densities of neutrons and protons separately, at 
LO = 19.4 + 0.5i MeV for ^^O and at w = 16.1 + 0.5i MeV 
for the K = mode of ^''Mg. These energies correspond 
to the peak energies of the GDR. In the transition density 
of ^^O, one can clearly see the isovector character. In the 
top-right panel of Fig. [U we show dpo/dr for ^^O, where 
the po{r) is the ground-state density profile. This corre- 
sponds to the transition density predicted by the Tassie 
model for an irrotational and incompressible fluid. As is 
discussed in the literatures [Si] , the GDR of light nuclei is 
approximately represented by this classical fluid model. 
On the other hand, the transition density of 2'*Mg shows 
more complicated features. We show, in the bottom-right 
panel of Fig. [SI a contour plot of the transition density 
in the Tassie model, 6p cx Vpo • VrFio(r) ~ dpo/dz. 
Compared to the RPA transition density, there is a de- 
viation, especially in the nodal structure in the interior 
region. This may suggest the importance of coupling of 
the dipole mode to the octupole and higher-multipole 
modes. 



C. Heavy nuclei 




E[MeV] 

FIG. 10: Calculated photoabsorption cross sections for ""Zr, 
^^"Sn, and ^'^^Pb. The calculation has been performed with 
Rbox ~ 15 fm, the SkM* parameter set, and 7 = 1 MeV. The 
experimental data (symbols) are taken from Refs. [stI. IssI. [39l| . 



splitting, however, this may be due to the spurious effect 
coming from the box discretization. Except for this split- 
ting, the results agree well with the experimental data. 
Table Hi] shows the calculated and the experimental val- 
ues of the GDR peak positions and the widths. Both 
values are extracted from the Lorentzian fit to the pho- 
toabsorption cross section. According to Ref. [Hj, the 
Lorentzian function of the form (I16p . 



a-abs(-B) = 



CTO 



(27) 



is used. The GDR peak positions are well reproduced 
within an error of 400 keV. The systematic deviation of 
the peak energies, seen in calculations for light nuclei, 
cannot be observed for these heavy systems. Therefore, 
we may conclude that the SkM* functional reproduces 
peak energies of the El resonances in heavy nuclei. 

The calculated widths are slightly smaller than the ex- 
periments, by 300 ~ 700 keV. Since we use 7 = 1 MeV in 
these calculations, this indicates the spreading width of 
= 1.3 ~ 1.7 MeV. The total damping width is about 
4 MeV for these nuclei. Thus, the spreading width is 
less than half of the total width. This is rather surpris- 
ing because, for heavy nuclei, the spreading width was 
suppos ed to be a major part of the total damping width 
|35l. l36l| . However, the fragmentation of the strength into 
non-collective Ip-lh states (Landau damping), which can 
be described by the RPA theory, is significant for these 
heavy systems. Thus, the spreading width of ^ 1.5 
MeV is able to reproduce a broadening of the experimen- 
tal strength distribution. This observation seems con- 
sistent with recent self-consistent RPA calculations for 
spherical nuclei 0] . 



For heavier nuclei, the calculation better agrees with 
experiments. As is well-known, in heavy spherical nuclei, 
a single energy-weighted Lorentzian curve can fit the ex- 
perimental data of the photoabsorption cross section very 
well [ill . We calculate photoabsorption cross sections in 
spherical nuclei ^°Zr, ^20Sjj^ a,nd 208p]3 within a box of 
radius -Rbox = 15 fm and compare them with experimen- 
tal data in Fig. (TUl The calculated GDR peak shows a 



IV. CONCLUSIONS 

We have presented an implementation of the finite am- 
plitude method (FAM) to make a fully self-consistent 
RPA calculation employing a realistic Skyrme energy 
functional. Although the RPA is a well established the- 
ory and have been widely applied, a fully self-consistent 
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TABLE II: Calculated and experimental [111 GDR peak en- 
ergies and widths for ^"Zr, ^^^Sn, and ^°®Pb. The FAM calcu- 
lation was performed with the same parameter set as Fig. 1101 





Calculation 


Experiment 


Nucleus 


Speak [MeVJ 


r [MeVJ 


-Epeak [MeVJ 


r [MeVJ 


«"Zr 


16.37 


3.85 


16.74 


4.16 


i2«Sn 


15.22 


4.52 


15.40 


4.89 


208 


13.26 


3.47 


13.63 


3.94 



calculation with a realistic functional is rather limited. 
This is mostly because of the complexity to construct the 
residual two-body kernels for realistic functionals. The 
FAM, which was proposed recently by the present au- 
thors, makes it possible to achieve a fully self-consistent 
RPA calculation without an explicit construction of the 
residual two-body kernels. Instead, the residual fields are 
evaluated by the finite difference, employing a computa- 
tional code for the static mean-field Hamiltonian alone 
with a minor modification. 

We implemented the FAM in the mixed representation 
where the orbital indices are used for the occupied or- 
bitals while the three-dimensional Cartesian grid points 
are used for the unoccupied orbitals. In this represen- 
tation, the linear response equation is a linear algebraic 
problem with a sparse and non-hermitian matrix. We 
have tested several solvers for the problem, and have 
found that the generalized conjugate residual method 
works stably to obtain the solution. We also examined 
the accuracy of the FAM employing a different param- 
eters of the finite difference, and have found that the 
method works robust, namely, accurate results can be 
obtained for a certain wide range of finite difference pa- 
rameter. 

We show results of the fully self-consistent RPA cal- 
culation for the electric dipole responses of several nuclei 
of spherical and deformed shapes. For light nuclei, we 
have found a systematic underestimation of the average 



excitation energy, irrespective of the energy functional 
employed. We further notice the underestimation of the 
strength above the giant dipole resonance. The discrep- 
ancy is significant only for light nuclei, the average exci- 
tation energies of heavy nuclei are reasonably described. 
For a deformed nucleus ^''Mg, the calculation shows a de- 
formation splitting of the giant dipole resonance, which 
well agrees with measurements. For spherical heavy nu- 
clei, we have found a substantial part of the width can 
be explained within the RPA, leaving less than half of 
the total width for a spreading width, in contrast to the 
previous studies which reported that most of the width is 
attributed to the spreading mechanism for heavy nuclei. 

We are currently performing a systematic calculation 
of the electric dipole responses for light to medium mass 
nuclei using the present method and Skyrme energy func- 
tionals. In a forthcoming paper, we will report a system- 
atic analysis of the properties of the giant dipole res- 
onance, including the average excitation energy, width, 
and the low-lying dipole mode around the threshold. 
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